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Given a video of M frames of size h X w. Background components of a 
video are the elements matrix which relative constant over M frames. In PCA 
(principal component analysis) method these elements are referred as 
“principal components”. In video processing, background subtraction means 
excision of background component from the video. PCA method is used to 
get the background component. This method transforms 3 dimensions video 
(hxw x M) into 2 dimensions one (N x M), where N is a linear array of 
size h x w. The principal components are the dominant eigenvectors which 
are the basis of an eigenspace. The limited memory block Krylov subspace 
optimization then is proposed to improve performance the computation. 
Background estimation is obtained as the projection each input image (the 
first frame at each sequence image) onto space expanded principal 
component. The procedure was run for the standard dataset namely SBI 
(Scene Background Initialization) dataset consisting of 8 videos with interval 


resolution [146 150, 352 240], total frame [258,500]. The performances are 
shown with 8 metrics, especially (in average for 8 videos) percentage of error 
pixels (0.24%), the percentage of clustered error pixels (0.21%), multiscale 
structural similarity index (0.88 form maximum 1), and running time (61.68 
seconds). 
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1. INTRODUCTION 

Background subtraction is an important step in many computer vision systems to detect moving 
objects [1]. It is commonly used in video surveillance applications to detect persons, vehicles, animals, etc., 
before operating more complex processes for intrusion detection, tracking, people counting, etc. The basic 
operation consists of separating the moving objects called “foreground” from the static information called 
“background”[2]. It consists in using a model of the scene background in order to detect foreground objects 
by differencing incoming frames with the model. Indeed, the first step in background subtraction is 
background estimation. 

We state the general problem of background estimation, also known as background initialization, 
bootstrapping, background reconstruction, initial background extraction, or background generation, as 
follows: Given a set of images of a scene taken at different times, in which the background is occluded by 
any number of foreground objects, the aim is to determine a model describing the scene background with no 
foreground objects [3]. 
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Background Estimation based on eigenspace was firstly proposed by Oliver et al [4]. Foreground 
don't appear in the same location in the M sample images and they are typically small, do not have a 
significant contribution to modeling background. Consequently, the portions of an image containing a 
foreground cannot be well-described by this eigenspace model (except in very unusual cases), whereas the 
static portions of the image can be accurately described as a sum of the various eigenbasis vectors. That is, 
the eigenspace provides a robust model of the probability distribution function of the background, but not for 
the foreground. 

The eigenspace model is formed by taking a sample of M images and computing both the mean and 
its covariance matrix. The eigenvectors of the covariance matrix are obtained by eigenvalue decomposition. 
In order to reduce the dimensionality of the space, in principal component analysis (PCA) [4]-[6] only 
k eigenvectors are kept, corresponding to the k largest eigenvalues. A CCTV (closed circuit television) can 
capture images in various resolutions, e.g., 176x120, or 352x240, or 704x480 pixels on the low to medium 
resolution range to 1280x720 or 1920x1080 pixels on the high resolution [7]. A CCTV with the lowest 
resolution, 176x120, size of the covariance matrix is 21120x21120. While computation eigenvalues and 
eigenvectors of large matrix take a long time [8]. 

The left singular vector at singular value decomposition of a normalized matrix which representation 
video was the principal component. Therefore, the singular value decomposition can yield the principal 
component without computing the covariance matrix. A subspace optimization technique to singular value 
decomposition significantly accelerates the classic simultaneous iteration method [9]. We propose Limited 
memory block Krylov subspace optimization for computing principal component. It will be used to construct 
background estimation. 


2. RESEARCH METHOD 

The following notation will be used: italics with subscripts to indicate vectors and matrix (Anw is a 
matrix of h rows and w columns), bold letters with subscripts for images (Bp, w is an image with height h and 
width w). Given an image I, of size h, w (height, width) it can be rearranged as a column vector xy,1, where 
N=h-w. 


1 column n” column 


v > 1“ column 
row > 
2" row > 
b n™ column 


Figure 1. Conversion of h x w image into N x 1 vector 


h? row > 


Given a video consists of M frames I7,I7,---,I@. These frames are arranged as a columns 


XN XN» Xa These vectors xy, 1, XA 1 °°» XN. are formed into a matrix X with height N and width M, 


1 2 M 
Xy, m = [XN XN XN | 
2.1. Principal Component via Limited Memory Krylov Subspace Optimization 
Eigenvectors of covariance matrix can be obtained by singular value decompositon. Suppose 


normalized matrix which representation video, Ay y = (x NM — LL" hm), where p is average image and J, m 
is row vector in which all elements are set to 1, singular value decomposition of matrix Ay y is 


Ay. = Uy, m=u,m Vun> 
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where Uy m is left singular matrix, Vm y is right singular matrix and Xm m is diagonal matrix in which 
diagonal elements are singular value. Multiplication Ay m with its transpose, Ajy,v, is given by: 


An mAMn = Uy mm, mVm, NVN mĒm,mUM,N> 
Cyn = Un.m@in MUM: (1) 


Equation 1 is eigenvalue decomposition of covariance matrix, Cy, where columns of Uy y are 
eigenvectors of Cy y. So, eigenvectors of covariance matriks are similar with left singular vector of Ay m. 
Eigenvectors of covariance matrix corresponding to the most significant eigenvalues will yield the principal 
components. So, principal component can be produced by compute left singular vector of matrix Ay m. X. 
Liu et all [9] proposed limited memory block krylov subspace optimization to compute dominant singular 
value decomposition. Its technique is proposed to significantly accelerate the simple subspace iteration 
method. The following we discuss how to obtain k largest singular value decomposition of matrix Ay y using 
limited memory block krylov subspace optimization. Starting from an initial point X € R“**, SSI (Simple 
Subspace iteration) computes the next iterate X*+” by the formula 


X@9D € orth(AA7X®); (2) 


where orth (M) denotes the set of orthonormal bases for the range space of M. As such, the iterates of SSI, 
with a possible exception for the initial guess, satisfy the orthogonality condition XTX = I. When k=1, the 
SSI method reduces to the well-known power method for computing the single largest eigenvalue and its 
eigenvector. In the SSI method, the orthonormalization step is indispensable (for example, see [10] for more 
details). 

For an unstructured matrix A, the computational costs of the matrix-block multiplication 
(i.e., AA’ X) and orthonormalization in (2) are O(MNk) and O(Mk?), respectively. In most applications, the 
approximating rank k is far less than the dimension m. Hence, the matrix-block multiplications of the type 
AA’ X constitute the dominate computational cost of SSI. Obviously, an acceleration will be achieved if one 
can reduce the number of iterations without having to incur extra matrix-block multiplications or other 
significant overhead. To achieve the goal of reducing the number of iterations, we propose to modify the 
basic SSI framework as follows. We replace the last iterate X® in the right-hand side of (2) by an 
“improved” intermediate iterate X so that 


X@®D € orth(AATX®), (3) 
where, for a chosen subspace S“ with a block Krylov subspace structure, 

X© := arg maxy-pnxkl|A7X||2, s.t.XTX =1, X ES, (4) 
Again, X € S® means all columns of X are from the subspace S. The selection of the subspace S which 
is constructed from a limited memory of the last a few iterates. Its choice is of course not unique. We first 
consider the subspace spanned by the current i-th iterate and the previous p iterates; i.e., 


sO = span (XO, XEN, <, X0); (5) 


where the memory length p > 0 will be specified in later. We collect the current and the other p saved iterate 
blocks in (5) into a matrix 


X = X® = [xO, X69,- , XCP] € RING; m 


where q = (p + 1)k is the total number of columns in x For notational simplicity, from here on we often 


choose to drop the superscript and subscript from quantities like xy whenever no confusion would arise. 
Also note that the collection matrix X is boldfaced to make it distinct from its blocks. Similarly, we define 


Y = YÊ = [yO, YOV,..., yE] € RM, o 
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which is also saved in memory. We emphasize that SSI already computes these blocks in Y and we only need 
to save them once computed. It is clearly that X € S{ if and only if X = XV for some V € R4**, and the 
subspace optimization problem (4) is equivalent to a generalized eigenvalue decomposition problem: 


max, <paxkll(ATX)V ||? s.t. V7(XTX)V = I. (8) 


However, numerical difficulty may arise in solving (8) as the matrix XTX can become numerically rank 
deficient. A more stable approach, which we will implement, is to find an orthonormal basis for 5, say, 


Q = Q” € orth(X®), 
and to express a matrix X € S® as X = QV for some V € R@**. Here we assume that X has a full rank and 


will later relax this assumption. We now convert the generalized eigenvalue problem (8) into an equivalent 
eigenvalue problem 


Maxycpaxk||RV||Z s.t. VTV =I, 9) 
where 
R =R; = A"Q). (10) 


Next we describe how to calculate the matrix product R in (10) from historical information without 
any additional computation involving the matrix A. Since Q € orth (X)and we assume that X has a full rank, 
there exists a nonsingular matrix C € R1%41 such that X = QC. Hence, Q = XC71, and R in (10) can be 
computed as 


R = ATQ = A'XC™ = YC, (11) 
Where Y = ATX is accessible from our limited memory. Once R is available, we can solve (9) by computing 
the k leading eigenvectors of the q x q matrix RTR. Let a solution to (9) be V. The matrix product in 
equation (3) can then be assembled as 


AATX® = ARV = AYC"'Y. (12) 


The remaining issue is how to efficiently and stably compute Q and R even when the matrix X is numerically 
rank deficient. We use the following procedure in our implementation. Noting that each block in X is 
individually orthonormal, we choose to keep the latest block X intact, and project the rest of the blocks 


onto the null space of X®", obtaining 
Py = PO = (1 — x0(x0)") [x69 = 2. xi], (13) 
Next we perform an orthonormalization of Py via the eigenvalue decomposition of its Gram matrix 
POP, =U, AU. (14) 


where Uy is orthogonal and A, is diagonal. It can be easily verified that the matrix 
i ; af ; 
Q =Q = x. Peden € orth PÊ, (15) 


provided that Ay is invertible. The above procedure can be stabilized by monitoring the numerical rank of 
Px, or specifically the eigenvalues on the diagonal of the matrixAy in (14). We choose to implement the 
following two-step stabilization scheme: 

Step 1 Delete the columns of Py whose Euclidean norms are below a threshold e, > 0. 

Step 2 Delete the eigenvalues in Ay, and corresponding columns in Uy, that are less than e€, > 0. 
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With a slight abuse of notation, we will continue to use Py, Uy and Ay to denote their stabilized versions, 
respectively, after possible deletions. Therefore, a stable construction of Q is still given by formula (15). 
After this stable orthonormalization, the corresponding R matrix can be generated as 


1 
RER = lr, Prost} (16) 
where before the stabilization procedure we had 


Py = PO = [yD ... yE] —~YO(xO)" [xen ws XE], (17) 


but some of the columns of Py may have been deleted corresponding to those deleted columns of Py during 
the stabilization steps. After the removal of numerical rank deficiency, the R matrix in (16) is well defined as 
is the Q matrix in (15). 

In summary, the algorithm performs eigenvalue decompositions on two small symmetric positive 
definite matrices: P{¥Py in (14) and RTR in (9). The sizes of the two matrices are pk and (p + 1)k, 
respectively, and frequently smaller due to deletions. Our computational experience indicates that in general 
p should be set to 2 or 3, or at most 4 but not greater. Consequently, when k is sufficiently smaller than N, it 
holds that (p + IKK M<N. 


2.1.1. Memory Length 

The memory length p, used for constructing the subspace S“ in (5), is a crucial parameter to the 
performance of our algorithm. The simplest way is to assign a constant integer value Pmax tO p at every 
iteration once the iteration counter i reaches Pmax; that is, at iteration i, 


p = min(i, Pmax). (18) 


In general a larger Pmax leads to a smaller number of iterations, but increasing Pmax also increases the 
computational costs per iteration. Our computational experiments indicate that usually a good balance is 
attained for Pmax E {2, 3, 4}. 

We have also found that an adaptive strategy on selecting p is useful to improving the performance 
of LMSVD. As the iterate sequence converges, the neighboring iterates tend to become more and more 
linearly dependent. Therefore, once judged appropriate it is beneficial to shrink the memory by deleting a 
block from the memory, reducing the size of later subspace optimization problems. Specifically, after Pmax 
iterations, we activate the following adaptive memory size strategy: 


ea eras ge (19) 


where [t] is the smallest integer greater than or equal to t, and Nc g) is the number of columns in R which 
can be smaller than(p + 1)k due to possible deletions done in the two stabilization steps. Combining (18) 
and (19), we reach our formula for selecting the memory length p at the i-th iteration: 


p=min fi |] =1, Pras (20) 


which is nonnegative. Generally, p initially increases to reach Pmax, then becomes non-increasing with a 
probability to decrease to a smaller value, even possibly to zero. Of course, when the memory length p 
becomes zero, our method reduces to the classic SSI. 


2.1.2. LMSVD Algorithm 
Based on the description above, we state our full Algorithm. For ease of reference, the algorithm 
will be referred to as LMSVD. 
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Algorithm LMSVD: Limited Memory Block Subspace Krylov Optimization for SVD 
1 Input Ay y and k. Initialize X = X® e€ R™**, Y = y% = ATX, andp = i = 0; 
2 while “not converged” do 

/* Block Subspace Optimization */ 

Compute Py by (13) and perform stabilization Step 1; 

Compute Py by (17) with the same column deletions as for Py; 

Compute the eigenvalue decomposition of PẸ Py in (14); 

Perform stabilization Step 2 to possibly shrink Xy and Uy; 

Compute R by (16) and eigenvalue decomposition of RTR; 

Let V solve (9), consisting of the k leading eigenvectors of RTR; 

Compute X = QV and Y = RV (which equals A7X); 

/* Simultaneous Iteration */ 

Compute XY € orth(AY®) and Y@Y) = ATX CD; 

Increment i, update p, X and Y, and continue. 

Output U, £, and V 


2.2. Background Estimation 

Once principal component (U) and the mean background (u) are computed, the input image (y;) 
with foreground objects was subtracted by the mean background. Defining an principal component matrix as 
U = [u1 U2," , Ug]. It follows that the coordinate (weight) in eigenspace of input image, y;, can be 
computed as follows 


w = (y; — u)" U, (21) 
when w is back projected onto the image space, a background estimation is created 
Yr = Uw! +u. (22) 


Noting that since the principal component matrix describes the general background appearances well 
however not the small moving objects, yg does not contain small objects. 


2.3. Metrics 

The metrics adopted to evaluate the accuracy of the estimated background models have been chosen 
among those used in the literature for background estimation [3]. Denoting with GT (Ground Truth) an image 
containing the true background and with CB (Computed Background) the estimated background image 
computed with one of the background initialization methods, the eight adopted metrics are: 

a. Average Gray-level Error (AGE): It is the average of the gray-level absolute difference between GT and 
CB images. Its values range in [0, L-1], where L is the maximum number of grey levels; the lower the 
AGE value, the better is the background estimate. 

b. Total number of Error Pixels (EPs): An error pixel is a pixel of CB whose value differs from the value 
of the corresponding pixel in GT by more than some threshold q (in the experiments the suggested value 
t=20 has been adopted). EPs assume values in [0; N], where N is the number of image pixels; the lower 
the EPs value, the better is the background estimate. 

c. Percentage of Error Pixels (pEPs): It is the ratio between the EPs and the number N of image pixels. Its 
values range in [0, 1]; the lower the pEPs value, the better is the background estimate. 

d. Total number of Clustered Error Pixels (CEPs): A clustered error pixel is defined as any error pixel 
whose 4-connected neighbors are also error pixels. CEPs values range in [0, N]; the lower the CEPs 
value, the better is the background estimate. 

e. Percentage of Clustered Error Pixels (pCEPs): It is the ratio between the CEPs and the number N of 
image pixels. Its values range in [0,1]; the lower the pCEPs value, the better is the background estimate. 

f. | Peak-Signal-to-Noise-Ratio (PSNR): It is defined as PNSR = 10: log; )((L — 1)*/MSE); where L is 
the maximum number of grey levels and MSE is the Mean Squared Error between GT and CB images. 
This frequently adopted metric assumes values in decibels (db); the higher the PSNR value, the better is 
the background estimate. 

g. MultiScale Structural Similarity Index (MS-SSIM): This is the metric defined in [11], that uses 
structural distortion as an estimate of the perceived visual distortion. It assumes values in [0; 1]; the 
higher the value of MS-SSIM, the better is the estimated background. 
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h. Color image Quality Measure (CQM): It is a recently proposed metric [12], based on a reversible 
transformation of the YUV color space and on the PSNR computed in the single YUV bands. It 
assumes values in db and the higher the CQM value, the better is the background estimate. 

While the last metric is defined only for color images, metrics 1 through 7 are expressly defined for 
gray-scale images. In the case of color images, they are generally applied to either the gray-scale converted 

image or the luminance component Y of a color space such as YcbCr . 


3. RESULTS AND ANALYSIS 

This section describes the results of the implementation of the proposed method. The result is 
background estimation using principal component. The metrics adopted to evaluate the accuracy of the 
estimated background models have been chosen among those used in the literature for background 
estimation. 


3.1. Background Estimation 

In this paper, the Scene Background Initialization (SBI) data set was chosen for the background 
estimation. The data set contains seven image sequences and corresponding ground truth (GT) backgrounds 
are given in Figure 1. In Table 1 we report, for each sequence, the name, the number of available frames, the 
subset of the frames adopted for testing, and the original resolution. The subsets have been selected in order 
to avoid the inclusion into the testing sequences of empty frames (frames not including foreground objects). 
The ground truths (GT) have been manually obtained by either choosing one of the sequence frames free of 
foreground objects (not included into the subsets of used frames) or by stitching together empty background 
regions from different sequence frames. Both the complete SBI dataset and the ground truth reference 
background images were made available through the SBMI 2015 website at http://sbmi2015.na.icar.cnr.it [3]. 


z., 
frame 0 frame 10 


ke 


CaVignal Foliage Foliage&People Snellen 


Hall&Monitor HighwayI HighwaylI 


Figure 1. Example frames from the seven sequences of the SBI dataset (first row) and coressponding GT 
(second row) 


Table 1. Information on Sequences Adopted for Evaluation 


Name Original frames Used frames Resolution 
Hall&Monitor 0-299 4-299 352X240 
Highwayl 0-439 0-439 320x240 
Highwayll 0-499 0-499 320x240 
CaVignal 0-257 0-257 200x136 
Foliage 0-399 6-399 200x148 
People&Foliage 0-349 0-340 320X240 
Snellen 0-333 0-320 146x150 


This algorithm is implemented on Matlab R2015a at personal computer with intel core i-5 processor 
and 4 GB RAM. the result of the implementation is the background estimation using the principal component 
via the LMSVD algorithm. The frame in the first row of Figure 1 is used as y; in equation 25, so we get the 
background estimation in Figure 2 using equation (26). The performance of the proposed method is evaluated 
both qualitative and quantitative. In Figure 2 we show the background images obtained by the proposed 
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methods on the SBI dataset, while in Table 2 we report accuracy results according to the metrics described 
in 2.3. 


Hall& 


Monitor Highwayl Highwayll CaVignal 
People& 
Foliage Foliage Snellen 
— n 


Figure 2. Background estimation, GT (Ground Truth) and CB (Computed Backgrounds were obtained by 
proposed method) 


Table 2. Accuracy Results of the Proposed Methods on the SBI Dataset 


Sequence AGE EPs pEPs CEPs pCEPs MSSSIM PSNR CQM 
Highwayll 3,24 273 0,00 1 0,00 0,98 33,24 32,92 
Highway! 3,71 2571 0,03 1400 0,04 0,95 29,82 29,82 
Foliage 3,87 1117 0,01 382 0,03 0,96 28,97 28,91 
Hall&Monitor 8,77 2849 0,10 2241 0,08 0,83 21,67 21,65 
CaVignal 12,14 4638 0,16 2778 0,09 0,95 24,32 23,74 
PeopleAndFoliage 34,93 47706 0,62 41913 0,54 0,74 15,24 14,76 
Snellen 43,55 15771 0,76 14416 0,69 0,75 14,25 14,16 
Average 15,74 10704 0,24 9019 0,21 0,88 23,93 23,71 


For sequence Hall&Monitor, a man walking straight down the corridor occupies the same image 
region for more than 65% of the sequence frames, while the briefcase is left on the small table for the last 
60% of sequence frames. The proposed method well handle the abandoned briefcase, but it includes a man 
who walks a bit obscenely. For both Highwayl and Highwayll sequences, the proposed method succeed in 
providing a faithful representation of the background model. This is due to the fact that, even though the 
highway is always fairly crowded by passing cars, the background is revealed for at least 50% of the entire 
bootstrap sequence length and no cars remain stationary during the sequence. 

For sequence CaVignal, the only man appearing in the sequence stands still on the left of the scene 
for the first 60% of sequence frames; then starts walking and rests on the right of the scene for the last 10% 
of sequence frames. The persistent clutter at the beginning of the scene leads of the proposed methods to 
include the man on the left into the estimated background. For sequence Foliage, even though moving leaves 
occupy most of the background area for most of the time, the proposed methods achieve a quite good 
representation of the scene background. 
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For sequence People&Foliage, the artificially added leaves and men occupy almost all the scene 
area in almost all the sequence frames. The proposed methods to include the contribution of leaves and men 
into the final background model. In sequence Snellen, the foreground leaves occupy almost all the scene area 
in almost all the sequence frames. The proposed methods achieve a quite good representation of the scene 
background. However, the scene background have low luminance than GT. 

In order to assess the challenge that each sequence poses for the proposed methods, we computed all 
metrics obtained by the proposed methods for each sequence, and ranked the sequences according to these 
metrics, as shown in Table 2. Here, Highwayl and HighwaylI sequences reveal as those that are best handled, 
while Snellen is the worst handled. Bearing in mind the kind of foreground objects included into the 
sequences, we can observe that their size is not a major burden; e.g., Foliage sequence is better handled than 
Hall&Monitor, even though the size of the foreground objects is much larger. As instance, CaVignal 
sequence is worse handled than Foliage, since it includes almost static foreground objects that are frequently 
misinterpreted as background. It can also be observed that the values of pEPs and MS-SSIM metrics 
perfectly vary according to the difficulty in handling the sequences; these two metrics confirm to be strongly 
indicative of the performance of background initialization methods. 

Conventional PCA can not be applied to get background estimation for the SBI dataset because the 
covariance matrix size of all sequences is too large as shown in Table 3. The hardware memory used is not 
sufficient to obtain the eigenvector of the covariance matrix (principal component). Table 3 shows the 
performance of the proposed method. The method succeed to produce the principal component for all 
sequence of images with time varies depending on total frames of sequences. 


Table 3. Performance of Proposed Method 


: Size of time 
Sequence Resolution covariance mains Total frames (proposed method) 
Hall&Monitor 352x240 84480x84480 295 frames 39,84 
HighwayI 320x240 76800x76800 440 frames 95,95 
Highwayll 320x240 7680076800 500 frames 223,83 
CaVignal 200x136 27200x27200 258 frames 6,15 
Foliage 200x148 29600x29600 394 frames 8,72 
People&Foliage 320X240 7680076800 341 frames 54,87 
Snellen 146x150 2190021900 321 frames 2,37 
Average 61,68 
Sequence Resolution ae of y Total frames ume 
covariance matrix (proposed method) 


4. CONCLUSION 

In this paper, we have presented a background estimation using the principal component. The 
principal component is obtained by compute the dominant singular value decomposition using the limited 
memory Krylov subspace optimization. The columns of left singular matrix of the dominant singular value 
decomposition is the principal component. Background estimation is obtained as the projection each input 
image (the first frame at each sequence image) onto space expanded principal component. 

The procedure was run for the standard dataset namely SBI (Scene Background Initialization) 
dataset consisting of 8 videos with interval resolution [146 150, 352 240], total frame [258,500]. The 
performances are shown with 8 metrics, especially (in average for 8 videos) percentage of error pixels 
(0.24%), the percentage of clustered error pixels (0.21%), multiscale structural similarity index (0.88 form 
maximum 1), and running time (61.68 seconds). 
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